package Core::Integration;
use utf8;
use strict;
use warnings;

our $VERSION = '1.00';
our @EXPORT_OK = qw(integrate extract_spectra calculate_mz);
use base qw(Exporter);

BEGIN {push @INC, '../'}

use JSON::XS;
use Core::Ascore qw(ascore);
use Core::IntConsts qw(get_masas_aa get_masas_mod);
use Parsers::Omssa;
use Parsers::Sequest;
use Parsers::Phenyx;

my %masas_aa = %{&get_masas_aa};
my %masas_mod = %{&get_masas_mod};

sub integrate{
	my @phenyx_files = @{$_[0]};
	my @sequest_files = @{$_[1]};
	my @omssa_files = @{$_[2]};
	my $database = $_[3];
	my $mgf_list = $_[4];
	my $report_file = $_[5];
	my $one_engine = $_[6];
	my $scan_method = $_[7];
	my $is_tmt_or_itraq = $_[8];
	my $offset = $_[9];
	my $precision = $_[10];
	my $experiment = $_[11];

	my %phenyx_info = Parsers::Phenyx::analyze(@phenyx_files);
	my %sequest_info = Parsers::Sequest::analyze(@sequest_files);
	my %omssa_info = Parsers::Omssa::analyze(@omssa_files);

	my %data;
	my @spectra_list;
	my @consensus_peptides;
	my %spectra;              #guarda nombres de espectros en keys, el valor no se usa
	my @todelete;

	foreach(keys %phenyx_info){$spectra{$_} = 1};
	foreach(keys %sequest_info){$spectra{$_} = 1};
	foreach(keys %omssa_info){$spectra{$_} = 1};

	print 'number of Sequest spectra:  ' . keys(%sequest_info) . ".\n";
	print 'number of Omssa spectra:  ' . keys(%omssa_info) . ".\n";
	print 'number of Phenyx spectra:  ' . keys(%phenyx_info) . ".\n";
	print 'number of TOTAL spectra:  ' . keys(%spectra) . ".\n";


	print "Joining search engines\n";
	foreach my $spectrum(keys %spectra){
		#print $spectrum;
		#Joining the info from different search engines and creating a consensus peptide
		if(defined $sequest_info{$spectrum}){
			$data{$spectrum}{'sequest'} = $sequest_info{$spectrum};
			$data{$spectrum}{'ms'} = $sequest_info{$spectrum}{'ms'};
		}
		if(defined $phenyx_info{$spectrum}){
			$data{$spectrum}{'phenyx'} = $phenyx_info{$spectrum};
			$data{$spectrum}{'ms'} = $phenyx_info{$spectrum}{'ms'};
		}
		if(defined $omssa_info{$spectrum}){
			$data{$spectrum}{'omssa'} = $omssa_info{$spectrum};
			$data{$spectrum}{'ms'} = $omssa_info{$spectrum}{'ms'};
		}

		my ($consensus, $consensus_mod, $consensus_z) = consensus($data{$spectrum});

		#Si solo uso un buscador, no borro ninguna key y copio los consensus de este buscador
		if ($one_engine){
			$data{$spectrum}{'consensus'} = $data{$spectrum}{$one_engine}{'unmodified_sequence'};
			$data{$spectrum}{'consensus_mod'} = $data{$spectrum}{$one_engine}{'modified_sequence'};
			$data{$spectrum}{'consensus_z'} = $data{$spectrum}{$one_engine}{'z'};
			push @consensus_peptides, $data{$spectrum}{$one_engine}{'unmodified_sequence'};
			push @spectra_list, $spectrum;
		}
		elsif ($consensus){
			$data{$spectrum}{'consensus'} = $consensus;
			$data{$spectrum}{'consensus_mod'} = $consensus_mod;
			$data{$spectrum}{'consensus_z'} = $consensus_z;
			push @consensus_peptides, $consensus;
			push @spectra_list, $spectrum;
		}
		#Deleting the non-consensus data if more than one search engine
		else {
			delete $data{$spectrum};
			delete $spectra{$spectrum};
			next;
		}

		#deleting spectra with double modification
		if($data{$spectrum}{'consensus_mod'} =~ /\,/){
			delete $data{$spectrum};
			delete $spectra{$spectrum};
			print $spectrum, "\t->\tDeleted because of double mod at N-terminus\n";
			next;
		}

		#Searching for the MS3->MS2 info
		if($data{$spectrum}{'ms'} == 3){
			my $file = get_ms2_file_from_ms3($spectrum, $scan_method);
			push @spectra_list, $file;
		}
		#Adding Experimental Conditions reference
		$data{$spectrum}{'experiment'} = $experiment;
	}

	print "number of spectra:  " . keys(%data) . ".\n";

	#Obtain the complete list of proteins containing the peptide
	print "Obtaining protein list\n";
	my $db_matches = search_fasta_db($database, @consensus_peptides);
	my %db_matches = %{$db_matches};
	foreach my $spectrum(keys %spectra){
		my $consensus = $data{$spectrum}{'consensus'};
		$data{$spectrum}{'consensus_proteins'} = $db_matches{$consensus}{'prots'};
		$data{$spectrum}{'consensus'} = $db_matches{$consensus}{'seq'};
	}

	#Obtain the spectra info (mass and intensity arrays and peptide experimental mass)
	print "Obtaining spectra\n";
	my %spectral_data = extract_spectra(\@spectra_list, $mgf_list);
	foreach my $spectrum(keys %data){
		if(defined $spectral_data{$spectrum}{'parent'}){
			$data{$spectrum}{'spectral_info'}{'intensity_array'} = $spectral_data{$spectrum}{'intensities'};
			$data{$spectrum}{'spectral_info'}{'mass_array'} = $spectral_data{$spectrum}{'masses'};
			$data{$spectrum}{'spectral_info'}{'parent_mass'} = $spectral_data{$spectrum}{'parent'};
		}
		else{
			print "Mgf info not found for the file $spectrum \n";
		}
		if ($data{$spectrum}{'ms'} == 3){
			my $file = get_ms2_file_from_ms3($spectrum, $scan_method);
			if(defined $spectral_data{$file}{'parent'}){
				$data{$spectrum}{'ms2_spectral_info'}{'intensity_array'} = $spectral_data{$file}{'intensities'};
				$data{$spectrum}{'ms2_spectral_info'}{'mass_array'} = $spectral_data{$file}{'masses'};
				$data{$spectrum}{'ms2_spectral_info'}{'parent_mass'} = $spectral_data{$file}{'parent'};
				$data{$spectrum}{'ms2_spectral_info'}{'ms2_file'} = $file;
			}
			else{
				print "MS2 info not found for the file $file \n";
			}
		}
	}

	#Obtaining the consensus deltamz
	foreach my $spectrum(keys %data){
		my $mz = calculate_mz($data{$spectrum}{'consensus_mod'},
							  $data{$spectrum}{'consensus_z'});
		$data{$spectrum}{'consensus_deltamz'} =  $mz - $data{$spectrum}{'spectral_info'}{'parent_mass'};
	}

	#Obtaining quantitative iTRAQ/TMT info
	my $myfunc;
	my $infokey;
	if($is_tmt_or_itraq eq "itraq"){
		$myfunc = \&analyze_itraq;
		$infokey = 'itraq_info';
	}
	elsif($is_tmt_or_itraq eq "tmt"){
		$myfunc = \&analyze_tmt;
		$infokey = 'tmt_info';
	}
	unless ($is_tmt_or_itraq eq "none"){
		my $info = '';
		my $low_offset = $offset - $precision;
		my $high_offset = $offset + $precision;
		foreach my $spectrum(keys %data){
			if ($data{$spectrum}{'ms'} == 2){
				if(defined $data{$spectrum}{'spectral_info'}{'intensity_array'}){
					$info = &{$myfunc}($data{$spectrum}{'spectral_info'}{'mass_array'},
									   $data{$spectrum}{'spectral_info'}{'intensity_array'},
									   $low_offset, $high_offset);
				}
			}
			elsif ($data{$spectrum}{'ms'} == 3){
				if(defined $data{$spectrum}{'ms2_spectral_info'}{'intensity_array'}){
					$info = &{$myfunc}($data{$spectrum}{'ms2_spectral_info'}{'mass_array'},
									   $data{$spectrum}{'ms2_spectral_info'}{'intensity_array'},
									   $low_offset, $high_offset);
				}
			}
			$data{$spectrum}{$infokey} = $info;
		}
	}

	#Obtaining Ascore info
	print "Calculating Ascore for phosphopeptides\n";
	foreach my $spectrum(keys %data){
		if ((defined $data{$spectrum}{'consensus_mod'})and
			(defined $data{$spectrum}{'spectral_info'}{'parent_mass'})
		   ){
			if($data{$spectrum}{'consensus_mod'} =~ /(\(21\)|\(23\))/){
				print $data{$spectrum}{'consensus_mod'}, "\n";
				my %ascore_data = ascore($data{$spectrum}{'consensus_mod'},
										 $data{$spectrum}{'spectral_info'}{'parent_mass'},
										 $data{$spectrum}{'consensus_z'},
										 $data{$spectrum}{'spectral_info'}{'mass_array'},
										 $data{$spectrum}{'spectral_info'}{'intensity_array'},
										 0.5, 2000, 200);
				$data{$spectrum}{'ascore_data'} = \%ascore_data;
			}
		}
	}
	my $json = JSON::XS->new;
	my $json_text = $json->pretty->encode(\%data);
	open my $SAL, ">", "$report_file";
	print $SAL $json_text;
	close $SAL;
	print "======================\n";
	print " Integration complete\n";
	print "======================\n";
	return 1;
}


sub get_ms2_file_from_ms3{
	my ($spectrum, $scan_method) = @_;
	my ($file, $scan1, $scan2) = split /\./, $spectrum;
	if ($scan_method == 0){                   #cid
		$scan1 = $scan1 - 1;
		$scan2 = $scan1;
	}
	elsif ($scan_method == 1){                #hcd
		$scan1 = $scan1 - 2;
		$scan2 = $scan1 + 1;
	}
	else{
		$scan1 = $scan1 - 4;                   #pqd
		$scan2 = $scan1 + 3;
	}
	return "$file.$scan1.$scan2";
}


sub consensus{
	my %data = %{$_[0]};
	#truco para la comparacion posterior: Doy valores por defecto diferentes
	my ($pept_sequest, $pept_omssa, $pept_phenyx) = ('0', '1', '2');
	my ($all_i_sequest, $all_i_omssa, $all_i_phenyx) = ('0', '1', '2');
	my $pept_cons;
	my $pept_cons_mod;
	my $z_cons;

	if(defined $data{'sequest'}{'unmodified_sequence'}){
		$pept_sequest = $data{'sequest'}{'unmodified_sequence'};
		($all_i_sequest = $pept_sequest) =~ s/L/I/;
	}
	if(defined $data{'omssa'}{'unmodified_sequence'}){
		$pept_omssa = $data{'omssa'}{'unmodified_sequence'};
		($all_i_omssa = $pept_omssa) =~ s/L/I/;
	}
	if(defined $data{'phenyx'}{'unmodified_sequence'}){
		$pept_phenyx = $data{'phenyx'}{'unmodified_sequence'};
		($all_i_phenyx = $pept_phenyx) =~ s/L/I/;
	}

	#no se puede poner el consensus_mode tan gratuitamente -> ascore!, peak assignation!
	if ($all_i_sequest eq $all_i_omssa){
		$pept_cons = $pept_sequest;
		$pept_cons_mod = $data{'sequest'}{'modified_sequence'};
		$z_cons = $data{'sequest'}{'z'};
	}
	elsif ($all_i_sequest eq $all_i_phenyx){
		$pept_cons = $pept_sequest;
		$pept_cons_mod = $data{'sequest'}{'modified_sequence'};
		$z_cons = $data{'sequest'}{'z'};
	}
	elsif ($all_i_omssa eq $all_i_phenyx){
		$pept_cons = $pept_omssa;
		$pept_cons_mod = $data{'phenyx'}{'modified_sequence'};
		$z_cons = $data{'phenyx'}{'z'};
	}

	return ($pept_cons, $pept_cons_mod, $z_cons);
}

sub search_fasta_db{
	my $fasta = shift;
	my @pepts = @_;

	my %matches;
	my %seq_match;
	my %prots;
	my %patterns;
	my $t0;

	$t0 = time();

	#elimino las secuencias repetidas para la busqueda
	my %nr_pepts;
	undef @nr_pepts{@pepts};
	my @nr_pepts = keys %nr_pepts;

	#calculo los patterns para que busque en database con I y L
	foreach (@nr_pepts){
		(my $pat = $_) =~ s/[IL]/[I|L]/g;
		$patterns{$_} = qr/$pat/;
	}

	$/ = "\n>";
	my $error = open my $FASTA, '<', $fasta;
	#Set cursor to start of file + 1 char (skip first ">")
	seek($FASTA, 1, 0);

	#search fasta
	while (<$FASTA>) {
		my ($title, $seq) = split(/\n/o, $_, 2);
		if ($title =~ /\|/){
			$title = (split /\|/, $title)[1];
		}
		else{
			$title =~ s/\>//;
		}

		$seq =~ s/\n//g;

		foreach (@nr_pepts){
			if ($seq =~ m/($patterns{$_})/){
				if (defined $matches{$_}){
					$matches{$_}{'prots'} = $matches{$_}{'prots'} . "|" . $title;
					push @{$seq_match{$_}}, $1;
				}
				else {
					$matches{$_}{'prots'} = $title;
					$seq_match{$_} = [$1];
				}
			}
		}
	}
	close $FASTA;

	#obtengo la sequencia con J si hay las dos alternativas I/L en la db
	while (my ($pept, $alt) = each %seq_match){
		if (@$alt > 1){
			my %unique;
			undef @unique{@$alt};
			my @unique_alt = keys %unique;
			if (@unique_alt > 1){
				my $first = shift @unique_alt;
				my @first = split(//, $first);
				for my $other (@unique_alt){
					my @other = split(//, $other);
					foreach (0 .. $#first){
						if ($other[$_] ne $first[$_]){
							$first[$_] = "J";
						}
					}
				}
				$matches{$pept}{'seq'} = join '', @first;
			}
			else{
				$matches{$pept}{'seq'} = shift @unique_alt;
			}
		}
		else{
			$matches{$pept}{'seq'} = shift @$alt;
		}
	}
	print 'time for database search:', time() - $t0, " secs\n";
	#init = 66 sec
	#after qr and $_ in 'for my $pept (@nr_pepts)' = 68 secs!
	return {%matches};
}


sub extract_spectra{
	#  Extract the parent mass and mass/intensity values
	#	@specs	->	list of spectra (file.scan1.scan2)
	#	@mgfs	->	list of mgf file inside we'll search for spectral info
	#	return	->	%specs{$spec}{'masses'} / {'intensities'} / {'parent'}
	my @specs = @{$_[0]};
	my @mgfs = @{$_[1]};
	my %specs;

	foreach (@mgfs){
		my %spec1 = parse_mgf($_, \@specs);
		foreach(keys %spec1){
			$specs{$_} = $spec1{$_};
		}
	}
	return %specs;
}


sub parse_mgf{
	my $mgf = shift;
	my @spec_list = @{$_[0]};
	my %specs_to_return;
	foreach(@spec_list){
		$specs_to_return{$_} = 1;
	}

	my $error = open my $MGF, '<', $mgf;
	my @lineas = do {local $/ = "END IONS\n\n"; <$MGF>};
	close $MGF;
	my %mgf_info;
	foreach my $linea(@lineas){
		$linea =~ /BEGIN IONS\nTITLE\=(.+)\nPEPMASS\=/;
		#iTRAQ_Nov09_EL_IMAC_Fr19,Scan:1015,MS:3,Rt:2220.280
		my $spec = (split /\,/, $1)[0];
		my $scan = (split /\,/, $1)[1];
		$scan =~ s/Scan\://;
		if($scan =~ s/\-/\./){
			$spec = "$spec.$scan";
		}
		else{
			$spec = "$spec.$scan.$scan";
		}
		if(defined $specs_to_return{$spec}){
			my @lineas_n = split /\n/, $linea;
			my $masas = '';
			my $intensidades = '';
			my $pept_mass = $lineas_n[2];
			$pept_mass =~ s/PEPMASS\=//;
			for(5..$#lineas_n){
				if($lineas_n[$_] =~ /END IONS/){
					last;
				};
				my ($mass, $int) = split / /, $lineas_n[$_];
				$masas = $masas . '|' . sprintf("%.4f", $mass);
				$intensidades = $intensidades . '|' . sprintf("%.2f", $int);
			}
			$masas =~ s/\|//;
			$intensidades =~ s/\|//;
			$mgf_info{$spec}{'masses'} = $masas;
			$mgf_info{$spec}{'intensities'} = $intensidades;
			$mgf_info{$spec}{'parent'} = $pept_mass;
		}
	}
	return %mgf_info;;
}


sub calculate_mz{
	my ($pep, $z) = @_;
	my $nt;
	my $mass = 18.010565;      #H2O

	while($pep){
		if($pep =~ s/^(\w)\((\w+),(\w+)\)//){
			$mass +=   $masas_aa{$1} + $masas_mod{$2} + $masas_mod{$3};
		}
		elsif($pep =~ s/^(\w)\((\w+)\)//){
			$mass +=   $masas_aa{$1} + $masas_mod{$2};
		}
		else{
			$pep =~ s/^(\w)//;
			$mass +=   $masas_aa{$1};
		}
	}
	return ($mass + $z * 1.007825) / $z;
}


sub analyze_itraq{
	my @masses = split /\|/, $_[0];
	my @intensities = split /\|/, $_[1];
	my $min = $_[2];
	my $max = $_[3];

	my ($int_114, $int_115, $int_116, $int_117) = (0, 0, 0, 0);

	for (1..$#masses){
		if($masses[$_] > 118){
			last;
		}
		#itraq ions -> 114.1112, 115.1083, 116.1116, 117.1150
		if(($masses[$_] >= (114.1112 + $min))and($masses[$_] < (114.1112 + $max))){
			$int_114 += $intensities[$_];
		}
		elsif(($masses[$_] >= (115.1083 + $min))and($masses[$_] < (115.1083 + $max))){
			$int_115 += $intensities[$_];
		}
		elsif(($masses[$_] >= (116.1116 + $min))and($masses[$_] < (116.1116 + $max))){
			$int_116 += $intensities[$_];
		}
		elsif(($masses[$_] >= (117.115 + $min))and($masses[$_] < (117.115 + $max))){
			$int_117 += $intensities[$_];
		}
	}
	return "$int_114|$int_115|$int_116|$int_117";
}


sub analyze_tmt{
	my @masses = split /\|/, $_[0];
	my @intensities = split /\|/, $_[1];
	my $min = $_[2];
	my $max = $_[3];

	my ($int_126, $int_127, $int_128,
		$int_129, $int_130, $int_131) = (0, 0, 0, 0, 0, 0);

	for (0..$#masses){
		if($masses[$_] > 131.6){
			last;
		}
		#tmt ions -> 126.1283 127.1316 128.1350 129.1383 130.1417 131.1387
		if(($masses[$_] >= (126.1283 + $min)) and ($masses[$_] < (126.1283 + $max))){
			$int_126 += $intensities[$_];
		}
		elsif(($masses[$_] >= (127.1316 + $min)) and ($masses[$_] < (127.1316 + $max))){
			$int_127 += $intensities[$_];
		}
		elsif(($masses[$_] >= (128.135 + $min)) and ($masses[$_] < (128.135 + $max))){
			$int_128 += $intensities[$_];
		}
		elsif(($masses[$_] >= (129.1383 + $min)) and ($masses[$_] < (129.1383 + $max))){
			$int_129 += $intensities[$_];
		}
		elsif(($masses[$_] >= (130.1417 + $min)) and ($masses[$_] < (130.1417 + $max))){
			$int_130 += $intensities[$_];
		}
		elsif(($masses[$_] >= (131.1387 + $min)) and ($masses[$_] < (131.1387 + $max))){
			$int_131 += $intensities[$_];
		}
	}
	return "$int_126|$int_127|$int_128|$int_129|$int_130|$int_131";
}


1;    #455
